Hydrogen Bonding and Vaporization Thermodynamics in Hexafluoroisopropanol‐Acetone and ‐Methanol Mixtures. A Joined Cluster Analysis and Molecular Dynamic Study

Abstract Binary mixtures of hexafluoroisopropanol with either methanol or acetone are analyzed via classical molecular dynamics simulations and quantum cluster equilibrium calculations. In particular, their populations and thermodynamic properties are investigated with the binary quantum cluster equilibrium method, using our in‐house code peacemaker 2.8, upgraded with temperature‐dependent parameters. A novel approach, where the final density from classical molecular dynamics, has been used to generate the necessary reference isobars. The hydrogen bond network in both type of mixtures at molar fraction of hexafluoroisopropanol of 0.2, 0.5, and 0.8 respectively is investigated via the molecular dynamics trajectories and the cluster results. In particular, the populations show that mixed clusters are preferred in both systems even at 0.2 molar fractions of hexafluoroisopropanol. Enthalpies and entropies of vaporization are calculated for the neat and mixed systems and found to be in good agreement with experimental values.


Introduction
The investigation of the mechanisms and interactions that regulate both neat fluids as well as their mixtures plays a fundamental role in the design of novel solvents. Over the last decades the interest in developing sustainable chemical processes has grown, and with it the need for new solvents that follow the fundamentals of green chemistry. [1] For example liquid-liquid -i. e. solvent -metal extraction offers a number of advantages compared to other techniques, [2] but the employed solvents are often toxic and expensive to dispose of. The last few years have seen an increasing demand of novel sustainable solvents for metal ions extraction. [3] Another approach is the use of the so-called deep eutectic solvents (DES), which are low-melting eutectics formed by the mixture of two or more substances whose eutectic point temperature is much lower than that of an ideal mixture. [4,5] With properties similar to those of ionic liquids, such as a low vapor pressure, low melting point, and high thermal stability, [6,7] they have come to be known as an economic and eco-friendly alternative for conventional organic solvents. A new kind of DESs, called Type V DESs, has recently been introduced. [8] It is formed by the combination of two non-ionic moieties which establish a strong hydrogen bond network, and has been proven to be more sustainable than traditional organic solvents. [8] In recent works, hydrophobic DES formed by menthol with different organic acids have been proved to present a strong hydrogen bond between the two components, stronger than the hydrogen bond network in the neat system. [9,10] Different hydrogen bond donors and acceptors have been investigated as components of hydrophobic or Type V DESs, [8] such as decanoic acid and lidocaine, [11] menthol with different natural acids, [12] and 1,1,1,3,3,3-hexafluoroisopropanol (HFIP) with betaine and L-carnitine. [13] HFIP in particular is an extraordinary solvent used in many different applications, including the activation of organic functionalities such as the intramolecular Schmidt reaction using a Lewis acid in HFIP, [14] the activation of carbonyl and epoxide substrates, [15,16] or the activation of hydrogen peroxide in the BaeyerÀ Villiger oxidation reaction. [17,18] It can also serve as a proton donor in dihydrogen bonding with different transition metal hydrides, as nonclassical hydrogen bond. [15,19] Its widespread use is due to a number of beneficial properties, such as its thermal stability, transparency to UV radiation, as well as its remarkable solvent properties, allowing it to dissolve a wide range of polymers, as well as most common polar and non-polar solutes. [20,21] In aqueous solutions, its acidity range is comparable to the one of formic acid. [15] Furthermore, its low boiling point facilitates its recovery via distillation. However computational studies on HFIP remain sparse. For instance HFIP has been investigated with molecular dynamics by some of the present authors to prove its catalytic effect on C,C coupling reactions on aromatic compounds with positive results. [16] In 2019, Deng et al. were the first to study HFIP-based DESs as the non-polar phase in liquid-liquid micro-extractions. [13] To be able to understand the mechanisms behind the formation of novel mixtures it is imperative to study their structure and in particular the hydrogen bond network which may exist in these liquids. From an experimental prospective, hydrogen bonds can be investigated via IR and NMR. [22,23] However, since their experimental detection and analysis can still be challenging and expensive, computational tools may be valuable for this kind of investigation. [24][25][26] We were able to calculate the activity coefficients of methanol in its binary mixtures with a set of small alcohols, which allowed some insight into how the chain length and branching modify the hydrogen bond network and consequently affect their behavior. [27] Since it is fundamental to study the hydrogen bond network in binary organic solvents, it is interesting to investigate the behavior of HFIP with both an hydrogen bond acceptor (acetone) and a compound able to work both as hydrogen bond acceptor and donor (methanol). The binary mixtures of HFIP with acetone or methanol have already been investigated experimentally in the past, but so far their liquid structure has never been investigated in detail. [16,28] In the current article, we investigate the inter-molecular interactions present in the binary mixtures of HFIP with either solvents. The study is based on a combined computational procedure of both classical molecular dynamics simulations (MD) and quantum chemical calculations. Furthermore, we use the binary quantum cluster equilibrium (bQCE) method to investigate the distribution of inter-molecular interaction motifs in the system. The method assigns populations to a set of clusters, enabling the weighting of properties according to the cluster distribution similar to Boltzmann weighting. In contrast to it, however, the bQCE method can weight clusters of different sizes and compositions and considers not only their electronic energies but also the particle volume and intercluster interactions. Boltzmann and bQCE weighting were compared in the past by our group in the field of the computational calculation of vibrational circular dichroism spectra of bulk phase, showing the advantages of our method over Boltzmann weighting. [29] For this purpose we use our inhouse code Peacemaker 2.8, which has proven to be a valuable tool to investigate the thermodynamic properties of both neat and mixed systems. [27,[30][31][32][33][34][35][36][37][38] The standard bQCE method involves optimization of two empirical parameters by fitting them to a set of experimental data such as boiling points or isobars; however, even for simple mixtures of organic solvents, these data are often unavailable in the literature. Hence we follow a different approach, namely, the isobars are instead obtained via MD simulations. In the following we will first outline and discuss the computational methods used, after which we will present and investigate the results we have obtained.

Computational details MD simulations
Initial configurations of the different systems were generated using the PACKMOL package (version 17.039). [39] Molecules were randomly placed in a cell with an initial cell vector of 40-50 Å, depending on the composition of the system. MD simulations were performed using the LAMMPS program package (version 11 Aug 2017). [40] OPLS-AA force field parameters were used for methanol (MeOH) and acetone, whereas for HFIP we adopted the force field developed by Fioroni et al., which was optimized to reproduce the experimental density. [41] A Lennard-Jones 6-12 potential was used to describe van der Waals (vdW) interactions. Lorentz-Berthelot mixing rules were used to describe non-bonded interactions. [42] A cutoff of 1 nm was applied for vdW and Coulombic interactions. In the first step the SHAKE algorithm [43] was used to constrain the bonds involving hydrogen atoms and followed by energy minimization. This process was repeated 3-5 times to let the systems mix correctly and eliminate energetic hot spots in order to stabilize the systems. Afterwards, the boxes were deformed to reach a preliminary and fixed cell volume, calculated to reflect a density of 0.8 g/cm 3 for each system. Following this, the systems were simulated for 1.5 ns in the NpT ensemble, using Nosé-Hover thermostat and barostat, to let the volume converge. The cell volumes over the last 0.5 ns were found to remain constant, as can be seen in Figure S2 in the SI. The system volume was then set to the average volume taken over the last 0.5 ns in the NpT ensemble. A further 1 ns of simulation time in the NVT ensemble was performed to further equilibrate the system. Finally, a production run of 20 ns was performed in the NVT ensemble using the same conditions as during equilibration. The time step was set to 0.5 fs for the pre-equilibration processes (shake, minimization, and deformation of the box) and increased to 1 fs afterwards. The simulations were analyzed with our in-house trajectory analysis code Travis. [44,45] The angular distribution functions (ADFs) were calculated using the cone correction included in Travis with a cut off of 250 pm as maximum distance between the reference and observed molecules. [46] Along with the Radial Distribution Functions (RDFs), the coordination numbers (CN) are calculated as well ðCN ¼ 1 bulk ∫g r ð Þr 2 drÞ, where g(r) is the RDF's intensity and 1 bulk is the bulk density of the system. Please note that, since the RDF is strictly dependent on the number of molecules that respect the given condition, the intensity of the systems cannot be compared. The peak's position, however, is not dependent on the number of molecules and they can be compared. Hydrogen bond lifetimes are analyzed using the dimer existence auto correlation function (DACF) implemented in the Travis code with the default curve fitting and approximations included in the code and described in literature. [47] For this analysis, the distance condition of 350 pm for the OÀ O distance and an angular condition of 135°-180°for the OÀ HÀ O angle were applied.

Cluster generation
The construction of a cluster set that is representative of the investigated system is a crucial step in the bQCE procedure and has been discussed in previous works. [27,32,38,48,49] As a first step, to find clusters, we performed a global minimum structure search for each cluster size and composition by running the genetic optimization algorithm OGOLEM [50,51] at the force field level of theory, using the AMBER 2016 program package [52] and the GAFF force field [53] implemented therein. The AMBER/OGOLEM combination is optimized to screen a great number of individual clusters. [27,32,38] The number of individuals per generation was varied between 80 and 240 in accordance to the cluster size to adjust for the increasing complexity. In total, a number of 2000-6000 individuals, i. e. clusters, were generated and evaluated for every possible composition up to a cluster size of six molecules. As the search for the global minimum structure by a genetic algorithm is performed on the classical force field level and the enormous configuration space poses a great challenge, the obtained structures are not necessarily identical to the global minimum on the quantum chemical level of theory. Instead, they can be understood as good candidate solutions to the global minimum and generally represent stable clusters that cover a range of enthalpically or entropically favored configurations. In addition, many of the obtained structures are expected to collapse towards the same geometry during quantum chemical optimization. Therefore, we select a number of ten clusters evenly distributed in the energetic range of the final generation for subsequent quantum chemical optimization for each cluster size. These clusters were then optimized at the BP86/def2-TVZP [54] level of theory with Turbomole (version 7.41) and an energy convergence threshold of 10 -9 . [55] The BP86 functional is proven to deliver good structural results and, if combined with a triple-ζ basis set, to reproduce with good accuracy the measured vibrational fundamentals. [56,57] The London dispersion energy was taken into account by applying the D3 dispersion correction. [58,59] Frequency calculations were performed for all clusters and those with imaginary frequencies were excluded. To avoid duplicate clusters in the cluster set, the structural similarity of all optimized clusters was quantified by their geometrical distance [60] d: (1) wherein I and I 0 are the principal moments of inertia of the clusters P and P 0 , respectively. Clusters P 0 with a geometrical distance of d P; P 0 ð Þ < 0:01 were removed from the cluster set. At the end of this procedure, 10 clusters of neat HFIP, 11 of neat acetone, 13 of neat methanol, 45 of the mixed HFIP-acetone system and 73 of HFIP-methanol were ready to be analyzed and to be included as inputs for the further bQCE calculations. Their interaction energy, size, and composition are tabulated in the supporting information. In the following, clusters will be given unique labels of the form hXsY-Z, where h stands for HFIP, s can be either acetone (a) or methanol (m), X is the number of monomers of HFIP, Y is the number of monomers of s, and Z is a label to differentiate clusters of the same composition.

The bQCE method
The theory of bQCE methods has been extensively detailed in several earlier works. [35,38,49,61] Through this method, we are able to calculate the cluster distribution and the thermodynamic properties of the system (for instance vaporization enthalpies) for a selected temperature range. Here, only a short overview of the bQCE method will be presented. As a first step, a system of noninteracting clusters in thermodynamic equilibrium is considered. Each cluster in that system is built up from either one (neat systems) or two (binary systems) monomers.
In thermodynamic equilibrium these clusters can transform into each other. We can write the equilibrium reaction between the clusters as where i P ð Þ and j P ð Þ are the number of monomers of each component C 1 and C 2 that form the cluster P. The system's total partition function Q tot at volume V and temperature T is given by where q tot P is the partition function of the single cluster P and N P f g is the full set of total cluster populations N P . This cluster partition function can be evaluated as product of partition functions corresponding to the different degrees of freedom: Here, q trans P is the translational, q rot P the rotational, and q vib P the vibrational partition function. They are calculated from standard equations for the particle in a box, rigid rotator, and harmonic oscillator, respectively. [61,62] The electronic partition function q elec P is calculated from the adiabatic binding energy D bind E elec P of the cluster. [63] Until here, we considered our system to consist of non-interacting clusters. However, in order to describe a liquid, not only the binding energy within a cluster but also the inter-cluster interactions must be considered. First, in order to take the volume into account that is taken up by the clusters themselves and is inaccessible to translation, an exclusion volume V ex must be subtracted from the phase volume V in q trans P . The exclusion volume is calculated as wherein b xv is the empirical exclusion volume parameter to correctly scale the cluster volume v P , which is sensitive to the choice of the volume method and the atomic radii used therein. In previous works, b xv was treated as temperature independent. [27,32,33,35,38] Here, we introduce a linear temperature dependence of where b xv is the exclusion volume expansion coefficient and b 0 xv is the base of the intercept. A similar approach was used in the past by Kelterer and coworkers. [64] Finally, the inter-cluster interactions must be taken into account. The electronic partition function q elec P is extended by a volume dependent mean-field-like correction term: where k B is the Boltzmann constant and the mean-field parameter a mf is an empirical parameter, that scales the strength of the intercluster interactions.
When performing a bQCE calculation, the empirical parameters b 0 xv , b xv , and a mf are chosen such that the deviation from a given reference property, such as densities and phase transition temperatures, is minimized. In earlier works, a simple grid sampling algorithm was used to optimize the empirical parameters. With the introduction of a third parameter b xv this method is no longer feasible. Here, the Differential Evolution algorithm [65] as implemented in the SciPy library [66] for Python 3.4 is interfaced with the Peacemaker 2.8 code to find the best solution. The script is available upon request.

Results and Discussion
Here, we will discuss the results of our investigation of the binary mixtures of HFIP with acetone and methanol, respectively. These were obtained by employing the bQCE method to a set of quantum chemically optimized clusters depicting different binding motifs, that are representative for the specific interactions in each system. Instead of employing experimental data, in this work, we use isobars obtained from a set of classical molecular dynamics simulations of the mixed systems that we conducted at different temperatures and mixture compositions as detailed in the computational details section.

Classical Molecular Dynamics Simulations
Classical molecular dynamics simulations of both mixed systems HFIP/acetone and HFIP/methanol were carried out at defined temperatures ranging from 298.15-338.15 K in intervals of 10 K. These simulations were repeated at different compositions with mole fractions of HFIP of 0.2, 0.5, and 0.8. The results are listed in Table 1 for HFIP/methanol (top) and HFIP/acetone (bottom). Both binary mixtures show similar behavior. Due to the high density of HFIP, the density increases with the mole fraction of HFIP but decreases with rising temperature. In addition, Table 1 lists the experimental densities of the HFIP/acetone mixture at 298.15 K, which were measured by Evans et al. using a single neck capillary tube pycnometer. [28] No experimental density could be found for the system HFIP/methanol. Our calculated densities for the system HFIP/acetone deviate by 2.5 %, 10.5 %, and 10.8 % from the experimentally measured values at mole fractions 0.2, 0.5, and 0.8, respectively. This deviation can have an impact on the thermodynamic properties of mixing calculated via the bQCE approach, however we demonstrate in the supporting information (SI) that a deviation in this range only slightly affects the cluster population and thermodynamical properties such as entropies and enthalpies of vaporization. Approximated thermodynamic properties of mixing are included in the SI for sake of completeness.
In order to characterize the hydrogen bonds present in the mixtures, additional analyses were carried out including Radial Distribution Functions (RDFs), Coordination Numbers (CNs), and Angle Distribution Functions (ADFs) of the simulated systems at 298.15 K averaged over 19.8 ns of production run. Figure 1 shows (from top to bottom) the RDF, the CN, and ADF of the OÀ H⋯O hydrogen bond between HFIP and acetone (left side) and between HFIP and HFIP (right side). Please note, HFIP is treated as the reference molecule. The observed molecule can be either acetone or another molecule of HFIP. The reference atom is the hydrogen of the hydroxy group, and the observed atom is the oxygen of the observed molecule. First, the inter-species hydrogen bond between HFIP and acetone will be considered. A sharp peak is present in the RDF at 161 pm at mole fractions of 0.2 and 0.5. The bond length decreases to 159 pm at a mole fraction of 0.8. An additional peak emerges around 350 pm at high concentrations of HFIP. The coordination number (CN) displayed in Figure 1 provides insight into the average number of molecules participating in a hydrogen bond. A single hydrogen bond is possible between the two molecules. At a distance of 250 pm, which is the location of the first minimum in the RDF and can be considered the maximum hydrogen bond length, the CN is 1.0 for the mole fraction of 0.2, meaning all the possible hydrogen bonds are filled by this interaction. With an increasing mole fraction of HFIP of 0.5 and 0.8, the CN decreases to 0.83 and to 0.26, respectively. In the RDF of the same-species hydrogen bond between two molecules of HFIP, a peak is present at 164 pm, which becomes sharper as the mole fraction of HFIP increases. At mole fractions of 0.5 and 0.8 an additional peak is visible at 326 pm, indicating an involved hydrogen bond network. From the CN plot, it is immediately clear that inter-species interactions are preferred over neat ones at mole fractions of 0.2 and 0.5, whereas same-species interactions become predominant at a mole fraction of 0.8. In particular, at a mole fraction of 0.2 same-species interactions are almost absent with a CN close to 0, but the CN increases to 0.17 and 0.73 at mole fractions of 0.5 and 0.8, respectively. This is due to a smaller number of acetone molecules able to be coordinated by the HFIP, forcing the same-species interaction to happen more frequently. The CNs for different binding motifs sum up neatly to 1.0, excluding the possibility of a significant presence of non-associated monomeric species in the system. Considering the CNs at 450 pm, which is the location of the second minimum in the RDFs, the inter-species interactions are preferred with a value of 1.3 against 0.9 of the H(HFIP)-O(HFIP) interaction at a mole fraction of 0.5. However at a mole fraction of 0.8 the situation is different, with the same-species interaction winning over the mixed ones with CNs of 2.2 and 0.6, respectively. A third peak is not visible at larger distances (RDFs up to 15 Å are included in the SI). This might indicate that circular configurations of neat HFIPs are present and maybe even preferred over chains at this mole fraction; however, this argument cannot exclude the presence of these linear formations. The bottom panels in Figure 1 show the ADFs of the hydrogen bond angle. Both the inter-species and same-species hydrogen bonds show similar behaviors and a clear preference for a linear arrangements.
In contrast to the HFIP/acetone mixture, both compounds in the HFIP/methanol can act as hydrogen bond donor or acceptor. Figure 2 shows the RDFs of the inter-species hydrogen bond length with HFIP acting as donor and acceptor, respectively. Looking at the (H)HFIP-(O)MeOH interaction first, respectively. The larger distance of 181 pm shows that the hydrogen bond donated by methanol is weaker than that donated by HFIP. Their position remains the same independent of the mixture's composition. The left and right side of Figure 3 show the RDF, CN, and ADF of the same-species hydrogen bonds shared between two molecules of HFIP and two molecules of methanol, respectively. The RDF of the (H)HFIP-(O)HFIP hydrogen bond shows a prominent peak at 161 pm and, similar to the same interaction in the mixture HFIP/acetone, as well as an additional peak at 324 pm at mole fractions of 0.5 and 0.8. Also similar to the previous mixture, the CN plot shows that only at mole fractions greater than 0.5 this interaction becomes predominant. In the RDF of the (H)MeOH-(O)MeOH interaction a peak is present at 180 pm. A second peak at 340 pm can be observed at mole fractions of 0.2 and 0.5. This peak is still weakly present at 0.8 around 380 pm. With an increasing mole fraction of HFIP, the same-species interaction of methanol with itself are concentrated in the first solvation shell, with a second solvation shell being present but less populated. The CN suggests that the MeOH-MeOH hydrogen bond is more significant in the equimolar mixture than the HFIP-HFIP hydrogen bond. As for the previous system, the CN values for the biding motifs sum up to 1. The CNs of the different interaction motifs sum up neatly to 1.0, which shows that non-associated monomeric species aren't present in significant amounts. Considering the CNs at 400 pm, which is the location of the second minimum in the RDFs shown in  oligomers cannot be ruled out, and there is proof in literature of chain configuration in pure methanol. [67] Overall, all the ADFs of this mixture at all mole fractions show that linear interactions are preferred.
The lifetimes τ of the hydrogen bonds in both mixtures were calculated as explained in the computational details and are presented in Table 2. A longer lifetime means a stronger bond, as it requires more time to break. In the top part of the table, the lifetimes of hydrogen bonds for the system HFIP/ Acetone are listed. The H(HFIP)-O(Ace) hydrogen bond has a lifetime of 9 ps, whereas the same-species H(HFIP)-(O)HFIP hydrogen bond has a significantly shorter lifetime of 3 ps. With an increasing mole fraction of HFIP, the lifetimes of the interspecies hydrogen bonds increase, meanwhile the same-species interaction remains stable. In the bottom part of Table 2, the hydrogen bond lifetimes for the system HFIP/MeOH are   [68][69][70] H(MeOH)-O(HFIP) is confirmed to be weaker than the other inter-species interaction.

Cluster analysis
Multiple clusters are optimized following the procedure described in the computational details section above to build the cluster sets which serve as input for the bQCE method. One of the most important features of the bQCE method is that it assigns populations to all clusters in the cluster set. The goal is to find the equilibrium distribution of all clusters so that they reproduce the reference isobars. The monomer-normalized population gives a measure of the importance of specific clusters and interaction motifs in the system. These populations are available for each temperature in the investigated temperature range.
In Figure 4 the most populated clusters in neat systems of pure HFIP, methanol, and acetone are shown for the temperature of 298.15 K. In the HFIP system, the tetramers h4-1 and h4-5 dominate the neat solvent, with populations of 0.57 and 0.37, respectively. In acetone the cyclic trimer a3-1 is the highest populated cluster with a population of 0.71. However, a significant amount of acetone molecules in the system exist in the form of non-associated monomers. The methanol system is dominated by ring formations, where the cyclic pentamers m5-1, m5-10, and hexamer m6-11 are the preferred geometries with values 0.37, 0.28, and 0.26, respectively. Figures 5 and 6 show visualizations of the most populated clusters in the binary mixtures HFIP/acetone and HFIP/methanol at 298.15 K, respectively. Since dozens of clusters have been   considered for both systems, here only the most populated ones are presented. Cartesian coordinates and visualizations of all cluster geometries are available in the supporting information. Focusing on Figure 5 first, at the mole fraction 0.2 the neat acetone trimer a3-1, which consists of three acetone molecules arranged in a triangular ring structure, is the most populated cluster with a value of 0.33, followed by the mixed hexamer h3a3-9 with 0.27. The same cluster, which can be described as aggregation of three alternating hydrogen bonded HFIP/ acetone dimers, is the most populated cluster in the equimolar mixture of HFIP and acetone with a population of 0.76. Please note that clusters with compositions that differ from the system's molar composition can be populated as the bQCE method will automatically find a cluster distribution that is consistent with the system's molar composition. In agreement with the CNs measured from MD simulations of the same system (Figure 1), there are no same-species bonds of HFIP in this cluster.
In contrast, at the higher mole fraction of 0.8, the neat HFIP tetramers are the most populated, which feature four HFIP molecules arranged in a quadratic ring of hydrogen bonds, each acting as both acceptor and donor. This indicates a mixture in which at lower concentrations of HFIP the interspecies hydrogen bond between HFIP and acetone is the most common interaction, in agreement with the CN plots in Figure 1. There is a good agreement between the CN at 400 pm and the RDF of the neat HFIP interactions and the most populated clusters at this mole fraction; already the classical simulation suggests the presence of neat HFIP rings at a mole fraction of 0.8. The distance between one H atom and the O atom on the molecule opposite to it in the tetramer is 327 pm, really close to the second RDF peak for the same-species interaction at 326 pm in Figure 1. At higher concentrations of HFIP, cooperative effects between multiple HFIP molecules become more important and the system is dominated by hydrogen bonded ring formations, while mixed interactions are still significant.
In Figure 6 the most populated clusters of the HFIP/ methanol mixture are displayed. At low concentration of HFIP, the highly coordinated hydrogen bond ring formations of the neat methanol pentamers dominate the system. But already at the low mole fraction of 0.2, the mixed h2m3-7 cluster, which contains hydrogen bonded HFIP and methanol molecules forming a ring in an alternating pattern, is highly populated. The possible presence of these kind of clusters was already suggested when discussing Figures 2 and 3. The significance of the mixed h2m3-7 cluster carries over to the equimolar mixture, where it is still the highest populated cluster and more populated than the stoichiometrically favored clusters h3m3-8 and h3m3-4, which are the second and third highest populated clusters. At high concentrations of HFIP, and in contrast to the HFIP/acetone mixture, the system is dominated by the mixed clusters h4m1-6, h4m1-7, and h4m1-8. This is in good agreement with the conclusions drawn from Figures 2 and 3.
Overall, there is a good agreement for both systems between the classical MD simulations and the populations calculated via the bQCE approach.
For now, the most significant clusters were considered only at room temperature. However, as the temperature changes other cluster formations might emerge and begin to dominate the bulk structure. The bQCE model provides such information and allows insight into the temperature dependence of the cluster equilibrium distribution in the selected temperature range of 298.15-338.15 K. Here, we will take a look at the populations of neat and mixed clusters and their evolution with temperature, which are shown in Figures 7 and 8 for different mole fractions of HFIP in the mixtures HFIP/acetone and HFIP/ methanol, where their population is significant (higher than 0.05). The top panel in Figure 7 shows the cluster populations in neat acetone. The trimeric a3-1 cluster dominates the system, but with rising temperature its population decreases in favor of the acetone monomer. In the center panel, cluster populations in the methanol system are presented. The pentamers m5-1 and m5-10 together with the hexamer m6-11 dominate the system until the boiling point is reached. Several research groups investigated neat methanol with the QCE approach. A recent work by Teh and coworkers, [71] which extensively investigates methanol cluster populations with both DFT and MP2 geometry optimization, states the most populated cluster is the octamer, which is not investigated in this article due to the computational time required, as explained in the computational details section. An older work by Kelterer and coworkers [64] finds instead that the most populated cluster size is the hexamer, followed by the pentamer, regardless that also clusters up to the octamer were investigated. The differences between those works and this article can be imputed to the geometry optimization and frequency calculation at different levels of theory, [71] or to a different procedure to generate the clusters. [64] However, there is a general agreement that multimolecular cyclic clusters are preferred until the boiling point, while in the gas phase the monomer dominates the population. In the bottom panel, the cluster populations in neat HFIP show that the tetramers h4-1 and h4-5 dominate the system until the boiling point. Figure 8 shows the cluster populations in the mixed systems at different mole fractions of HFIP. At low concentrations of HFIP, the neat acetone trimer is the preferred cluster at 298.15 K, whereas the monomer population increases significantly with rising temperature. This is due to the breaking of the enthalpically favored inter-species hydrogen bonds (calculated interaction energy of À 53.6 kJ/mol) and the rise of entropically favored small clusters such as the acetone dimer (calculated interaction energy of À 26.8 kJ/mol) and the non- associated acetone monomer. Pure HFIP clusters are not populated even at higher temperature. Instead, HFIP is bound in mixed clusters, of which the hexamer h3a3-9 is the most populated. In the equimolar HFIP/acetone mixture, the hexamer h3a3-9 is highly populated at 298.15 K, but with rising temperature decreases in population in favor of the acetone monomer and smaller mixed clusters dominated by HFIP. At high concentrations of HFIP, neat acetone clusters are almost absent. Meanwhile, the tetrameric neat HFIP clusters h4-1 and h4-5 are significantly populated, which is in agreement with earlier observations of the emergence of cooperativity effects at high HFIP concentrations. With rising temperature, the population of mixed clusters increases.
The right column of Figure 8 shows the cluster populations in the HFIP/methanol. In contrast to the other mixture, the methanol-dominated mixed cluster h2m3-7 is more populated at low concentrations of HFIP and the populations are nearly constant over temperature. The equimolar mixture is again composed mainly by mixed clusters, in particular the methanoldominated h2m3-7, whereas neat clusters are absent. At high concentrations of HFIP, the HFIP-dominant mixed cluster h4m1-6 dominates the system, but with rising temperature its population decreases in favor of pure HFIP clusters. Neat HFIP clusters are significantly populated, but less than in the HFIP/ acetone mixture.
For both systems, the formation of close inter-species interactions is favored over interactions of the same species. This behavior is more pronounced for the HFIP/methanol mixture.
To complete the quantum cluster equilibrium analysis, the interaction energies, the hydrogen bond lengths, and the hydrogen bond angles of the dimers, are presented in Table 3, wherein for mixed dimers the first named species is the hydrogen bond donor and the last named species is the hydrogen bond acceptor. It is apparent that the interaction energies of the mixed dimers, where HFIP acts as the donor, are significantly stronger than that of any neat dimer. Both the bond lengths and interaction energies of the isolated dimers indicate the importance of inter-species interactions over samespecies interactions if HFIP is the hydrogen bond donor, in line with the MD results presented before, where there are exclusively H(HFIP)-O(Ace) and H(HFIP)-O(MeOH) interactions at mole fractions of HFIP of 0.2 and 0.5.

Thermodynamic properties of neat and mixed systems
Through calculating the system's total partition function based on the equilibrium distribution of a set of representative clusters, the bQCE method grants access to the absolute thermodynamic functions such as the Gibbs energy G, enthalpy H, and entropy S at any investigated temperature. Using the bQCE method, we can thus calculate properties such as the enthalpy and entropy of vaporization D vap H and D vap S as simple difference of the liquid phase and a gas phase reference. Already in earlier works we were able to establish our procedure of using a so-called QCE 0 calculation, wherein a mf is set to 0 removing all inter-cluster interactions, as gas phase reference. [27,32,72] Table 4 compares the calculated vaporization enthalpies and entropies of the neat systems to their experimental reference values at room temperature and the boiling point. Overall, good to excellent agreement with the experimental reference is achieved for all systems. This is true both at room temperature and at boiling points of the solvents. In all cases, D vap H is slightly overestimated, possibly indicating an over-stabilization of the liquid phase. The largest deviation is observed for methanol at its boiling point. A likely explanation is the experimentally observed aggregation of methanol molecules to small clusters in the gas phase, [73] which is not properly sampled by the QCE 0 calculation, as the monomers are populated with 99 %. This is in agreement with our already published calculated D vap H of 39.33 kJ/mol of methanol at room temperature presented in a previous work, calculated at the same level of theory. [27] Here, we obtain a slightly different value due to the changes in methodology.
Enthalpies and entropies of vaporization at 298.15 K were calculated for the mixed systems as well, see Table 5. For HFIP/ acetone the calculated vaporization enthalpies are higher than those of the neat components, with the highest value calculated for a mole fraction of 0.5. In contrast, the entropies of vaporization are always higher than that of pure acetone, but lower than that of pure HFIP. In HFIP/methanol at a mole fraction of 0.2 the calculated enthalpy of vaporization is lower than those of the neat components. At higher mole fractions, the enthalpy of vaporization is higher than those in either of the neat components, similar to the HFIP/acetone system. The same behavior is observed for the entropy of vaporization.

Conclusions
The mixtures of HFIP with methanol and acetone were investigated. First, molecular dynamics simulations were performed to get the isobars of the systems for the temperatures 298.15, 308.15, 318.15, 328.15, and 338.15 K at mole fractions of HFIP of 0.2, 0.5, and 0.8, respectively. Then the structure of the hydrogen bond from the simulations were analyzed and discussed. Molecular dynamics simulations and hydrogen bond evaluation from the dimer clusters calculated at DFT level are in good agreement. Together with the population analysis and the evaluation of the most populated clusters, we are presenting two strongly interacting mixtures, with some hints that the HFIP/methanol system presents a stronger hydrogen bond framework. Clusters of up to six molecules were optimized at DFT level for both the mixtures and the neat systems. The hydrogen bonds of the dimers were analyzed. The calculated clusters, from simulated isobars and experimental data, where available, were used as input for the binary quantum cluster equilibrium method to get the thermodynamical properties. The enthalpies and entropies of vaporization of the neat systems are in good agreement with previous theoretical results and with experimental values. In addition, enthalpies and entropies of vaporization were calculated for the mixed systems.